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Recent experiments have shown that the striking structure formation in dewetting films of evaporating col- 
loidal nanoparticle suspensions occurs in an ultrathin 'postcursor' layer that is left behind by a mesoscopic 
dewetting front. Various phase change and transport processes occur in the postcursor layer, that may lead to 
nanoparticle deposits in the form of labyrinthine, network or strongly branched 'finger' structures. We develop 
a versatile dynamical density functional theory to model this system which captures all these structures and may 
be employed to investigate the influence of evaporation/condensation, nanoparticle transport and solute trans- 
port in a differentiated way. We highlight, in particular, the influence of the subtle interplay of decomposition in 
the layer and contact line motion on the observed particle-induced transverse instability of the dewetting front. 



How surface patterns and structures evolve over time is of 
great interest for a wide range of scientific fields. Striking 
examples include river network patterns (T), the growth of 
rocks around geothermal springs [2] evaporation-caused cof- 
fee stain patterns 0, and the patterns in the distribution of 
living organisms (4). Many such structures are generated by 
the interaction of fluid motion over the surface and deposition 
and/or abrasion of material. A particular process of high re- 
cent interest that concerns us in this Letter is the formation 
of structures during the (evaporative) dewetting of nanoparti- 
cle suspensions on solid substrates (5H2). The patterning is 
generic to a wide class of dewetting evaporating suspensions 
and solutions 18UT21 and depends crucially on the interplay of 
several competing phase-change and transport processes. 

The rapidly expanding study of such systems currently re- 
ceives strong impetus from research in two distinct areas. 
On the one hand, studies from the last couple of decades of 
the dynamics of dewetting of surfaces by non- volatile liquids 
lfl~3l ITU have been extended to investigate the interplay be- 
tween (de) wetting and evaporation of volatile liquids (T5l[T6ll . 
On the other hand, there is interest in the non-equilibrium 
thermodynamics and rheology of the respective phase and 
flow behaviour of bulk suspensions and solutions - see e.g. 
Ref. ifTTl and references therein. A thin film of pure non- 
volatile liquid that is deposited upon a smooth substrate (e.g., 
a polystyrene film having a thickness of a few tens of nanome- 
ters, deposited on silicon oxide [13]) may rupture due to ef- 
fective molecular interactions between the film surface and the 
solid substrate. The rupture mechanism can be (i) via a sur- 
face instability (often called spinodal dewetting) that occurs 
spontaneously and results in patterns of a certain characteristic 
wavelength, or (ii) via nucleation at randomly distributed de- 
fects (HDDS). The resulting holes then grow to form a polyg- 
onal network of liquid rims that may subsequently decay into 
drops. All stages of this process are intensively studied: the 
rupture mechanisms (13] [331 [T8J, the hole growth (20), the 
morphologies and evolution of the resulting patterns [21 ], and 
the stability of receding liquid fronts [ 22, 23 ]. For reviews of 
this body of work, see Ref. [24]. 

The dewetting processes of solutions and suspensions are 
more involved than those of a pure liquid because they in- 
volve several interdependent dynamical processes: transport 
of solute or colloids, transport of the solvent and evapora- 



tion/condensation of the solvent. As a consequence, one 
must distinguish between 'normal' convective dewetting and 
evaporative dewetting. Experimental studies performed with 
volatile solutions/suspensions of polymers (U |9), macro- 
molecules (T0llT2l and colloids a few nanometers in size (re- 
ferred to as 'nanoparticles') (5lf71 [251 describe a variety of 
richly structured deposits of the solutes. One may observe 
labyrinthine and polygonal network structures similar to the 
structures observed following 'classical' dewetting. As the 
solvent evaporates, the solute remains dried onto the sub- 
strate and therefore 'conserves' the transient dewetting pattern 
(SlfTOl. However, the solute is not just a passive tracer: it may 
influence the thresholds and the rates of the initial film rup- 
ture processes. Most importantly, it may also destabilise the 
straight dewetting fronts and trigger the creation of strongly 
ramified structures - a process observed in many different sys- 
tems (7U9ll26l. 

For instance, the dewetting of films of a suspension of thiol- 
coated gold nanoparticles in toluene may result in the depo- 
sition of the nanoparticles in branched finger patterns. The 
precise properties of these depend on the strength of the attrac- 
tion between the colloidal particles Q. Employing contrast- 
enhanced video microscopy to study the dynamics of the sys- 
tem, initially the receding of a mesoscopic dewetting front 
(equivalent to the receding three phase contact line) is ob- 
served, leaving behind an unstable ultrathin 'postcursor' film, 
having a thickness similar to the diameter of the nanoparticles. 
Subsequently, the postcursor film breaks, forming a pattern of 
holes that themselves grow in an unstable manner, resulting 
in an array of branched structures. Note that the mesoscopic 
front may also be unstable - an effect that does not concern us 
here. 

Theories for modelling such processes are rather limited at the 
present time. Hydrodynamical models of dewetting by thin 
films based on a long- wave approximation (24J |27) provide a 
mesoscale description of the liquid film, and can account for 
evaporation of volatile liquids (28l and for the presence of a 
solute (29l . However, they are not able to describe the pro- 
cesses occurring in the postcursor film as they do not account 
for the interactions between the solute particles and between 
the solute particles and the solvent. Alternatively, to describe 
these processes one may employ two-dimensional (2d) kinetic 
Monte-Carlo (KMC) lattice models that focus solely on the 
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FIG. 1: (Color online) (a) Phase diagram of the pure liquid in the 
plane spanned by the density pi and temperature T. (b) Liquid and 
nanoparticle densities pi and p n at coexistence as a function of the 
average nanoparticle density p n , for eu/ksT = 1.25, eni/ksT = 

0.6, Snn — 0. 



dynamics of the solute diffusion and the solvent evaporation 
l30H32l . The neglect of convective solvent transport can be 
justified based on estimates comparing its influence with the 
one of evaporative solvent transport l32l . However, so far the 
2d KMC models have not incorporated diffusive solvent trans- 
port that might be important in the postcursor layer. 
In this Letter we present an alternative description for the 
structure formation in the postcursor film that does not have 
the limitations of the previous approaches. We develop a 2d 
dynamical density functional theory (DDFT) 133H361 to de- 
scribe the coupled dynamics of the density fields of the liq- 
uid pi(Y,t) and the nanoparticles p n (r, t). In this approach, 
diffusive liquid transport can be incorporated in a straight- 
forward manner enabling us to go beyond previous 2d KMC 
studies and examine its influence. To construct the DDFT 
model, we (i) develop via coarse-graining an approximation 
for the underlying free energy functional of the system and 
(ii) form equations governing the dynamics of the two den- 
sity fields. These are able to account for the non-conserved 
and conserved aspects of the dynamics, i.e., phase change and 
diffusive transport processes, respectively. 
In order to compare results with the established KMC lattice 
model [30-32], we start from the lattice Hamiltonian to de- 
velop a mean-field (Bragg- Williams) approximation for the 
free energy functional (35j [37). Expressing the interaction 
terms (sums over neighbouring lattice sites) as gradient opera- 
tors l35lL the following semi-grand l38l free energy functional 
is obtained: 



F\PhPn] = J 



dr 



+£nl(^p n ) ■ (Vpi) ~ ppl 



Pi)] 



where f(pup n ) = k B T[pi\npi + (1 - p/)ln(l 
k B T[p n \np n + (1 - p n ) ln(l - p n )] - 2supf - 2s nn p n - 
^niPnPh includes entropic contributions and various inter- 
action terms - the parameters eij, where z, j — n,l, are the 
energies for having neighbouring pairs of lattice sites occu- 
pied by species i and j, respectively, T is the temperature, k B 
is Boltzmann's constant and we have set the lattice spacing 
(7 = 1. Note that F in Eq. ([I]) can also be obtained by making 
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FIG. 2: (Color online) Results for spinodal evaporative dewetting 
of a nanoparticle suspension, (a) and (b) are typical nanoparticle 
density profiles, for times t/ti — 6 and 20, and (d) and (e) are the 
corresponding liquid density profiles, for Mf = 0. The domain size 
is 200cr x 200cr. (c) and (f) give the corresponding time evolution 
of the mean liquid density (pi) and the structure factor (p n (k) 2 ) for 
various values of Mf. The remaining parameters are £u/kBT = 



1.25, £ni/k B T = 
and (p n ) — 0.3. 



0.6, s ri 



0, a = 0.4M^ nc cr 4 , p/k B T = -3.4, 



a gradient expansion of the (non-local) free energy functional 
of a continuous system [37]. The rate of evaporation of the 
liquid from the substrate is determined by the chemical poten- 
tial p in the reservoir (i.e. in the vapour above the substrate). 
When the temperature is sufficiently low, and the chemical 
potential p = p coe ^ we observe coexistence between a thick 
(high density) and a thin (low density) liquid film. In Fig. 1(a) 
we display the limit of linear stability (spinodal) and the equi- 
librium coexistence curve (binodal) for the pure liquid (i.e. 
with p n = 0). To indicate the influence of the solute we plot 
in Fig. 1(b) the densities pi and p n at coexistence for a fixed 
temperature k B F jeu = 0.8 as a function of the average con- 
centration p n = \ (p n +p n ), where p n and p h n are the densities 
of the nanoparticles in the coexisting a and b phases. For fur- 
ther details concerning the phase diagram and its topology see 
Ref. |39l. 

At equilibrium, the derivative p n = 5F[p ni pi]/5p n is a con- 
stant, corresponding to the chemical potential of the nanopar- 
ticles. However, when the system is out of equilibrium, p n 
may vary along the substrate. We assume that the thermo- 
dynamic force V// n drives the dynamics of the nanoparticles 
and that the nanoparticle current is j = — M n p n V p n , where 
M n (pi) is a mobility coefficient. This expression for the cur- 
rent, together with the continuity equation, yields the time 
evolution equation for the nanoparticle density profile: 



dp n 

dt 



M nPn V 



5F[p n ,pi] 



Spn 



(2) 



This equation may also be obtained by assuming over-damped 
stochastic equations of motion for the nanoparticles I33ll34l . 
To model the fact that nanoparticles do not diffuse over the 
dry substrate (when pi is small) we set the mobility M n (pi) to 
switch at pi = 0.5 (smoothly) from zero for the dry substrate 
(low pi) to a for the wet substrate (high pi). Note that our 
results are not sensitive to the precise form of M n (pi). 
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FIG. 3: (Color online) Density profiles for evaporative dewetting in 
the nucleation regime, (a)-(c) are the nanoparticle density profiles at 
times t/ti = 20, 80 and 4000, (d) and (e) are the liquid profiles for 
t/ti = 20 and 80, for Mf = 2 and p/k B T = -3.33; remaining 
parameters are as in Fig. [2] In (f) we plot the average density of the 
liquid on the substrate (pi), as a function of time for this case and 
the case Mf = 0. The system was initialised with the (discretised) 
density profiles: pi(x, y, t = 0) = 0.9 + 0.05%, p n (x, y,t = 0) = 
0.3 + 0.27%, where % is a random number uniformly distributed on 
the interval [—1,1]. It is due to this random noise that the holes are 
nucleated in some places and not in others. 
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FIG. 4: (Color online) Density profiles from the evolution of an 
unstable dewetting front, (a)-(c) are the nanoparticle density pro- 
files at times t/U = 2000, 20000 and 40000, (d) and (e) are the 
liquid profiles for t/U = 2000 and 20000, for a = 0.2Mf c a 4 , 
p/ksT = —3.28 and domain size 800a x 800cr; remaining param- 
eters are as in Fig. [2] The initial density profiles are the same as those 
in Fig. [5] except for y < we set pi = p n = 10 -10 , to create an ini- 
tially straight dewetting front at y = 0. In (f) we plot the mean finger 
number (/) as a function of the mobility coefficient for nanoparticle 
diffusion a, for the case when Mf — 0. 



For the liquid, the density may change either by evapora- 
tion/condensation from/to the substrate (non-conserved dy- 
namics) or may diffuse over the substrate (conserved dynam- 
ics). The latter dynamics is treated in a manner analogous 
to that for the nanoparticles. For the non-conserved dynam- 
ics, we assume a standard form (40) , i.e., that the change 
of the density over time is proportional to — (/i sur f — p) = 
— SF\p n , pi]/Spi, where /i S urf(r,£) is the local chemical po- 
tential of the liquid on the substrate. Combining these two 
contributions, we obtain the time evolution equation of the 




FIG. 5: (Color online) Nanoparticle density profiles from the evolu- 
tion of an unstable dewetting front, for the time t/ti = 20000, (a) 
for Mf = and (b) Mf = 5, for a = 0.5Mf c a 4 , p/k B T = -3.28 
and domain size is 800cr x 800cr; remaining parameters and initial 
density profiles are as in Fig. [4] In (c) we show the dependence of 
the mean finger number (/) on the mobility coefficient for liquid 
diffusion, Mf. 



liquid density profile: 



dpi 
dt 



Mf Pl V 



SF{p n ,pi] 



5 Pi 



M? c 



SF[p n ,pi] 
dpi 



(3) 



We assume that the two mobility coefficients Mf and M z nc are 
constants. In what follows we set = 1 and Mf 10 = 1. 
Note that in the low density limit, when p n —> and p\ — » 0, 
the conserved part in both Eq. ^ and Eq. ^ corresponds to 
Fickian diffusion. 

Before discussing results from our theory, we first make a cou- 
ple of comments about the status of the theory. Firstly, we note 
that Eq. ([]]) constitutes a simple 'zeroth-order' mean-field ap- 
proximation for the free energy of the system and omits (for 
example) terms such as ln(l — p n — pi) which describe the 
excluded area correlations between the liquid and the nano- 
particles. Second, due to the fact that we derive the theory 
from the (already) coarse-grained lattice Hamiltonian rather 
than by integrating over degrees of freedom (coarse-graining) 
in the full DDFT theory for the three-dimensional liquid film, 
the theory can not be regarded as a 'fully' microscopic theory. 
The reason that we have modelled the system using this sim- 
ple theory is because our interest is the basic question of what 
physics drives the behaviour displayed in the experiments. In 
this work we choose to start from the lattice theory, in order to 
compare with the KMC. However, as future work we plan to 
go beyond the lattice theory. Our goal here is not to construct 
a model describing every detail of these systems; instead we 
seek to examine what physics is involved in determining the 
observed pattern formation. Nonetheless, this DDFT does al- 
low us to investigate the time evolution of the postcursor film 
of an evaporating nanoparticle suspension under fewer restric- 
tions than the KMC model. To discuss the importance of liq- 
uid transport within the layer we compare results obtained for 
different liquid mobilities (Mf > 0) and results without liquid 
transport (i.e., setting Mf = 0). We focus on three examples: 
(i) spinodal dewetting (see Fig. [2]), (ii) dewetting via nucle- 
ation of holes in an initially flat film (see Fig. [3J, and (iii) the 
unstable receding of an evaporative dewetting front, which ex- 
hibits branched fingering (see Figs. [4] and [5}. 
Fig. [2] shows snapshots from a purely evaporative spinodal 
dewetting process. Panel (c) gives the evolution of mean liq- 
uid density with time, (pf) = A~ x J drpi(r,t), where A is 
the area of the substrate, and panel (f) gives the nanoparticle 
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structure factor S(k) = (p n (k) 2 ), where p n (k) is the Fourier 
transform of p n ( r )- For small times, the unstable film de- 
velops a typical spinodal labyrinthine pattern with a typical 
wavelength 2tt / /c max (note that the symmetry breaking in the 
density profiles in all our calculations is due to the addition of 
a small amplitude random noise to the density profiles at time 
t = and no noise is added at later times). The nanoparticles 
concentrate where the remaining liquid is situated. However, 
they are 'slow' in their reaction: when pi already takes val- 
ues in the range 0.08 to 0.83, the nanoparticle concentration 
has only deviated by about 25% from its mean value. The 
film thins rapidly forming many small holes. The competi- 
tion for space results in a fine-meshed network of nanoparti- 
cle deposits with a much higher concentration of particles at 
the network nodes - an effect that can not be seen within the 
KMC model. Because the liquid wets the nanoparticles, some 
liquid always remains on the substrate. Accounting for sol- 
vent diffusion, the rate of the dewetting process is increased 
[see Fig. |2jc)] and leads to a more strongly modulated final 
pattern - i.e. the peaks in S(k) are higher for Mf > than for 
Mf = Qf)]. 

Fig. [3] shows snapshots from a dewetting process triggered 
by nucleation events. The holes nucleate at several arbitrary 
places due to the random initial noise in the density profiles, 
and grow to form a random polygonal network of rims of 
highly concentrated solution. On a very long time scale the 
network coarsens into an array of drops. The influence of liq- 
uid transport can be seen in the final panel of Fig. [3] where 
we display a plot of the average density of the liquid on the 
substrate as a function of (log) time. We see that the liquid is 
able to evaporate from the substrate faster and that for times 
t/ti ~ 10 2 — 10 4 , the total amount of liquid on the substrate is 
less when Mf = 2, than when Mf = 0. However, over very 
long times, the drops on the substrate slowly move and 'eat- 
up' the network pattern. This process is faster with solvent 
diffusion, leading to a faster increase of (pi) at long times. 
Since the liquid wets the nanoparticles, some liquid also re- 
condenses back onto the substrate. 

The final example in Figs. [4] and [5] is the evolution of the fin- 
gering instability for a receding dewetting front. The fingering 
instability is caused by a build up of the nanoparticles at the 
receding front, which collects the nanoparticles due to their 
attraction to the liquid. In Fig.|4ja) we see that at early times 
the initially straight front shows a rather short-wave instabil- 
ity; about 20 short 'fingers' can be seen. However, the fin- 
ger pattern coarsens rapidly to a stationary pattern containing 
only about half the initial number of fingers. Intriguingly, the 
mean finger number remains constant although at the mov- 
ing contact line new branches are created and old branches 
merge continuously. The occurrence of this phenomenon in 



the present continuum model (DDFT) is similar to results of 
the KMC [32], and proves that jamming of discrete particles 
is not a necessary mechanism for causing the instability. In 
Fig.|4jf) we show how the average number of fingers per unit 
length, (/), varies as a function of a, the mobility coefficient 
of the nanoparticles on the wet substrate, for the case when 
Mf = 0. We see that as a is decreased, the number of fin- 
gers increases. This increase in (/) occurs because when the 
mobility of the nanoparticles is decreased the front 'collects' 
more particles (less of them diffuse further from the front). 
The resulting region of high concentration solution at the front 
may be 'dynamically unstable': As the front velocity depends 
non-linearly on the amount of particles collected, any fluctua- 
tion along the front may trigger a transverse instability. 
Fig. [5] shows that the finger number (/) depends non- 
monotonically on the mobility coefficient for liquid diffusion, 
Mf. Although the overall trend is an increase of (/) with in- 
creasing Mf, there exists an intermediate region (5 < Mf < 
7) where (/) slightly decreases. The overall trend results from 
an increase in front velocity (due to the increase in Mf) at 
fixed particle diffusivity. However, we currently have no ex- 
planation for the intermediate slight decrease. 
Note also that in all cases the instability may be strongly am- 
plified when the particle interactions favour the clustering of 
the nanoparticles (this is not the case for the parameters used 
here), where the higher concentration at the receding front 
leads to a local demixing of nanoparticles and liquid, that it- 
self enforces the deposition of a highly branched finger pat- 
tern. In this 'demixing regime' the instability is determined 
by the dynamics and the energetics of the system whereas for 
the case studied in Fig.|4]it mainly depends on the dynamics. 
We note finally, that the fingering process can be seen as a 
self-optimisation of the front motion, so that the average front 
velocity is kept constant by expelling particles into the fin- 
gers. A similar effect exists for dewetting polymer films (22j, 
where surplus liquid is expelled from the growing moving rim 
which collects the dewetted polymer. However, front instabil- 
ities found for dewetting polymers only result in fingers with- 
out side-branches flTTl or fields of droplets left behind l22ll . 
In this Letter we have developed a versatile DDFT that is ca- 
pable of describing the pattern formation observed in evapo- 
rating dewetting thin films of suspensions. Since our DDFT 
takes account of all the basic solvent and solute transport and 
phase change processes in a consistent manner, we believe 
that it will be the basis for many successful future studies of 
the behaviour of suspensions and solutions at interfaces. 
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